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SUMMARY 


Several iterative algorithms based on multigrid methods are introduced for solving linear 
Fredholm integral equations of the second kind. Automatic programs based on these algorithms 
are introduced using Simpson’s rule and the piecewise Gaussian rule for the numerical 
integration. 


INTRODUCTION 


Several multigrid iterative methods based on the Nystrom method are applied for the fast 
solution of the large dense systems of equations that arise from the discretization of Fredholm 
integral equations of the second kind. We will consider the linear Fredholm integral equation of 
the second kind, 

Ax(s) - f k(s,t)x(t)dt — y(s), seD (1) 

J D 

with D a bounded close domain, and yE X where X is the underlying Banach space. Necessary 
assumptions are 

(i) k(s , t ) is such that the associated integral operator K is compact from X into X 

(ii) A is not an eigenvalue of I\ and A ^ 0 


The Nystrom method for solving (1) uses some type of numerical integration to obtain the 
approximating equation 

Xxi(s) ~^2a j (s)x t (t j ) = y(s), sED (2) 

i = i 

the nodes t u t 2l are in D, and x t (t) = x(t). The weights a^s) can be defined in a variety of 

ways, depending on the smoothness and form of the kernel function. If k(s,t) and x(t) are 
reasonably smooth, usually aj(s) = wjk(s,tj ), where 
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( 3 ) 


f f{t)dt « J2 w jf( t i) 

D j = 1 

is a numerical integration formula. Let the numerical integration operator Ki be defined by 


7l[ 


Kix(s) = s£D 
i=i 

( 4 ) 

Using (2) and (4), (1) approximated by the linear system 


71, 

= y(t,) 

j=i 

( 5 ) 

We will denote (1) and (5) symbolically as 


(A - K)x = y 

and 

( 6 ) 

(A - Ki)xi = y 

( 7 ) 


respectively. Our discussion is based on the convergence of a sequence of approximations to the 
unique solution of (1). 


In finding numerical solutions for equations (1), the system (5) is too large to be solved 
directly. The purpose of this paper is to consider some iterative variants of (4). The basic 
assumptions needed in our algorithms are given in section 2. In section 3, linear iterative 
algorithms are given based on Simpson’s rule and piecewise Gaussian quadrature rule for the 
numerical integraion formulae. And in the section 4, we include numerical examples. 


BASIC ASSUMPTIONS 


The methods will be defined and discussed using the abstract formulation of Anselone [1] and 
Atkinson [3], [4] for families of collectively compact operators. 

Let X/, l = 0, 1, 2...,be finite-dimensional subspaces of the Banach space X and let 
P/, / = 0, 1, 2, ..., be a bounded projection operator from X onto X/. We need the following 
assumptions for {X;} and {P,} 

(Al) X 0 cX l C .... C X,... C X 
(A2) lim ||/ - Pif || = 0 for all f 6 X 

I—+OO 
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The sequence {X/} is thought as being associated with a sequence of decreasing meshsizes {h;} 
with lim hi — 0. Corresponding with this sequence {/q},we approximate K by a sequence of 

operators {A';}, Ki : X — * X. In multigrid iteration, the subscript l is called ’Tevel”.The 
hypotheses on {Ki : l > 1} and K are as follows. 

(A3) K and K[ , / > 1 are linear operators on the Banach space X into X. 

(A4) Kix — > Kx as n — ■* oo, for all x € X. 

(A4) {Ki} is a collectively compact family of operators. 

The following is a consequence of the assumptions (A3) - (A5): 

Lemma 1 Assume (AS) - (A5). Then with n defined as in (3) 

(i) K is compact 

(ii) | \(K — Ki)K\\ and ||(A" — Ki)Ki\\ converge to zero as n oo 
(Hi) If at = sup sup II (AT — K m )K n \\ , then lim a/ = 0 

m>( n>l , - >0 ° 

Proof. See Atkinson [4], 

Lemma 2 //(A — A') -1 exists, then 

(A — K[)~ l exists for sufficiently large say N( A), and is uniformly bounded by c 2 (A) and 

\\x-x,\\ <c 2 (X)\\Kx- K,x\\, 1>N( A) 

where xi = (A — Ki)~ l y 
Proof. See Atkinson [4], 

This shows xi — » x and gives a rate of convergence. 


LINEAR ITERATIVE METHODS 


Multigrid Methods 


Assume that x/, o denotes a approximate solution of (7) with residual 

di = yi-(X~ Ki)xi fi ( 8 ) 

Then improve on the accuracy by writing 

%i, i = xi i0 + hi (9) 

where the correction Si satisfies the residual correction equation 
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(A - K,) 6 , = d, 


( 10 ) 


In general, the correction term 61 will be small, and it is unnecessary to solve the residual 
correction equation (10) exactly. Thus we may write 


61 = Bidt (11) 

where Bi denotes a bounded linear operator approximating (A — K{)~ 1 . By (??) and (9) together 
with (11), we obtain 


xi , i = [A - Bi ( A - Ki)]xi fi + Biy t (12) 

as the new approximate solution to (7). The equation (11) can be represented well by means of 
coarser grid functions 


(A — AVi)6/_i = i (13) 

where is chosen reasonably and depends linearly on d/. If r : Xi — * Xi_ i is the restriction 
mapping, then 

di-i = rdi (14) 

Having defined d/_iby (14), £;_i is obtained using (11) at level /— 1. Having obtained <5/_xwhich 
is defined only on the coarse grid level, we need to interpolate this coarse-grid function by 

St = p6i- 1 (15) 

where p describes the prolongation of a coarse grid function to a fine grid function. 

We note here that the choice of the prolongation p in (15) must be small enough to satisfy 

\\I-pr\\<ChJ (16) 

where the consistency order r depends on the discretization, (e.g. on the order of the quadrature 
formula). For the restriction operator r, we will consider both trivial injection and Nystrom type 
restriction. 

Our automatic algorithm is based on the following multigrid iteration which is given as a 
recursive procedure. 


Multigrid iteration for solving (A — Ki)x{ = y 
Procedure Multigrid (/, x/,y) 

if / = 0 then 

solve xq = (A — Ki)~ l y 
otherwise 


(17) 
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XI = \[I<ixi + y] 
di = (/ - Ki)xi - y 
di_i = rd\ 

repeat the Procedure Multigrid with (/ - 
x? ew = Xl “ p^(-i 

We now give some basic results of the multigrid algorithm (17) that are used in our automatic 
algorithm. 


Let Cfc be the contraction number of the multigrid iteration employed at level k 


-i+i 


~Xk < (k 


x J k - x k 


(18) 


Then it is known that {(*} are uniformly bounded by some C < 1- 


Let 


C := max (k (^) 

where l is the maximum level in (17). The relative discretization error, the difference between x k 
and Xk- 1 , is often estimated by 


\\pxk-i - Zfc|| < Ci h\ 

for 1 < k < / 

where p is a prolongation operator and t is the consistency order. 


( 20 ) 


Theorem 3 Assume (20) and 


C 2 C < 1 


( 21 ) 


with 

h^y 

hk 

then the i th iteration of the multigrid procedure (17) at level k results in x k and satisfies the 
error estimate 


Co := max 

i <k<i 


Xk — %k\\ < CzC\h k T 
for 0 < k < l 


where 


C 3 = 


c 

1 - C 2 c* 


( 22 ) 


(23) 
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Proof. See Hackbush [11]. 


Theorem 4 Assume the validity of (22) and suppose < - then the i th iteration of the 
multigrid procedure (17) at level k results in Xk satisfies the error estimate 

(24) 

(25) 

Proof. See Hackbush [11]. 


where 


Xk ^ C 4 


c 4 = 


(2 T — 1)C 

1 - C 2 c* 


Automatic Algorithms 


The automatic algorithm £t in (18) is used to estimate the iteration error. Then together 
with the discretization error the global error in the solution is estimated. Often £* is estimated by 


(k = 


„j + 1 

x k 


4 


X 


i-i 


Then 


x k - 1 


a 


1 - (k 


J + 1 


~ XI 


(26) 


(27) 


is used to estimate the iteration error. Thus at any level, a minimum of two iteration is required 
to estimate the iteration error. However, (24) together with (25) can be used to estimate £ using 


c 4 = 


iteration error 


discretization error 
and it will enable us to estimate (27) with only one iteration. 


(28) 


Our first algorithm is based on Simpson’s rule with double the node points as the level 
increases, i.e. dimension of the linear system at a level l is 2 ,+1 + 1. In this case we have C 2 = 16 
in (21). Thus by the condition (21), if £ < ^ the estimates in (22) holds with i=l, i.e. only one 
multigrid iteration per level. The result is computational savings. As the level increases the 
amount of computation increases, so that there is a significant time savings in performing only one 
iteration as the dimension of the linear system being solved becomes larger. Moreover £ k in (18) 
goes to zero as the level k increases, which means that after a certain level k, £t becomes so small 
that the iteration error becomes much less significant than the discretization error, hence more 
accurate estimation of it is not needed. Thus one iteration is sufficient at this stage. 
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The second algorithm is based on the piecewise Gaussian quadrature rule for the numerical 
integration scheme. We adapt the iteration error estimation scheme discussed earlier. 

For simplicity we use hi = for l — 1, 2, ... This means that we reduce the length of each 
subinterval by half as the level increases. Suppose at some level /, we have a partition 


Ql = {a = q 0 <q r < .... < q mi = b} 


(29) 


with 


q, = a + i*h t for i = 0, 1, 2, ..., m/ 
and m t = 2‘ :=number of subintervals, for / = 0, 1 , 2, .... 


Then 


rb ™t v 

J a - , 


i=l i= 1 


where 


/ = f(h) 

JO • i 


(30) 


(31) 


is the Gaussian quadrature rule on [0,1] with p node points. 

Unlike Simpson’s rule, we do not have nested node points. In the following algorithm, both 
restriction and prolongation are done with Nystrom type interpolation. 

Procedure Multigrid with piecewise Gaussian (/, x;,y) (32) 

if / = 0 then 

solve x 0 = (A — A'o) - 1 2/ 
otherwise 

x, = \[l<ixi + y] 

di = (A — Ki)x( — y — K\xi - Kixi 
di - 1 = r(K t x t - K t xi ) 

repeat the Procedure Multigrid with (l- 1 ,^_ 1 ,<? ; _i) 
x?' w = x t - p8i_ x 


Nystrom type interpolations as in the procedure (32) are costly. Each interpolation involves 
O(rii) multiplications at each level. However this can be improved as suggested in our conclusion 
later. 

The following theorem which is due to Atkinson-Potra [7] gives the theoretical iterative rate of 
convergence for piecewise Gaussian quadrature with Nystrom type interpolation. We will assume 
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that the kernel k(s,t) belongs to the class G(a,”f). This means that the kernel k(s,t) has the 
following properties: 

(Gl) Define 

= {(s,t)\a<$<i <b} 

^2 = {(s ,0 | a < t< s < 6 } 

Then there are functions k, £ C a {'&i),i = 1,2 
with 

k(s,t) =^(3,1), (s,t) € # 1 , t?8 

k{s,t)=k a(M). (M )€*2 

(G 2 ) If 7 > 0, then fc(s,t) € C 7 ([a, 6 ] x [a, 6 ]). If 7 = -1, then the kernel k(s, t ) may have a 
discontinuity of the first kind along the line t s 

Theorem 5 Assume that k(s,t ) € G{a, 7 ). Then solve the Nystrom equation 

N 

X({s ) = Y w iH s i + ^( 5 ) 

i = 1 

using piecewise Gaussian quadrature rule with p node points in subintervals by first 
obtainning xfiti), xi(tiv) as a solution of the linear system 

N 

xi(U) = Y Wjk(ti,tj)xi(tj ) + y{U) ( 34 ) 

j=l 

then using (33) as an iterpolation formula gives an error estimate 

||x = X,|| = O(V) (35) 

where w — min {a, 2 p ,7 + 2 }. 


Proof. See Atkinson-Potra [7] for the case p=r+l. 

Finally to determine i, the needed number of iteration at any level l, use (24) and (25) with 
t == 2 p, hence C 2 = 2 2p . 


Automatic Implementation 


Our automatic implementation is divided into two stages based on the results from the 
iteration method. In stage 1, (A - K m )x m = y is solved directly, and then an attempt is made to 
solve (A - K,)x, = y for l > m, iteratively. If the rate of convergence is sufficiently rapid then 
the stage 2 is entered. Otherwise m is replaced by l and the stage 1 is repeated. In stage 2, the 
value of m will serve as the coarsest grid level in the multigrid procedure (17) and solve 
( X _ K,)x, = y iteratively until termination of the algorithm. The iteration procedure attempts to 
use the minimum number of iterates such that once the iterative solutions satisfy a certain criteria 
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we will try to estimate the rate of convergence asymptotically, which enables the estimation of the 
rate of convergence with only one iteration per level. As shown in our numerical examples, this 
scheme results in computational savings at finer grid levels. 


The initial guess for an iteration of the higher level is the interpolation of the solution of the 
preceding level which may have been obtained either directly or iteratively. The error ||x — x m || 
and ||ar — xj|| in stages 1 and 2, respectively, are monitored continuously, regardless of whether the 
iteration method is being used or not. Thus the multigrid iteration may not have been invoked 
successfully before the attainment of an answer within the desired error tolerance. 


In order to estimate the global error in the current solution, we need to monitor the 
discretization error and the iteration error. For the iteration error estimation, (27) is used with 
estimated £ in place of £*,. In stage 1, a test is made to determine whether the speed of 
convergence is sufficient to enter stage 2. If 

C < [Ratio] 1 / 2 (36) 

then the speed of convergence is adequate for stage 2. This requirment will usually insure that 
only two iterates are needed to be calculated in stage 2 at any given level. The number Ratio is 
the theoretical rate at which the error in xi should decrease when l is increased to the next level. 
In our case, since we are doubling the node points as the level increases, Ratio = with r = 4 
for Simpson’s rule and r = 2p for p points piecewise Gaussian quadrature in each subinterval. 


For the discretization error estimation, we compute the rate at which the error is decreasing 
for the current level. For each computed level /, 


NumDE := ||x,-x/_i|| (37) 

and let DenDE be the previous value of NumDE, if any. Then the rate is computed using 


DE;= NumiyE 

DenDE 

Using this value of DE, we estimate the error x — x/, 


(38) 


Error := 


DE 1 
.1 - DE. 


NumDE 


(39) 


which is a standard error estimate for sequences which are converging geometrically with a rate 
DE. Having estimated Error as in (39), we use the final test 


Error < t (40) 

with 6 a desired error tolerance supplied by the user. 


To ensure that only needed accuracy in Xi is computed, we want to test 

iteration error < quadrature error (41) 
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This is done by 


,( 2 ) _ mi 




l i 


< 


w 
, c 


) (r^f) IN 2 ’ - x! °’l 


The test (42) is obtained by using (41) and the approximations 





(42) 


(43) 


X — X 


( 2 ) 



(44) 


If the test (42) is not satisfied, then the new iterate is calculated, and (42) is tested again. 
Once an iterate is acceptable according to (42), we check for accuracy in the most recently 
computed iterate using (39) and (40). 


NUMERICAL EXAMPLES 


The integral equation 


rb 

;(s) — A I k(s,t)x(t)dt = y(s), a < s < b 
J a 


is solved with the kernel 


(45) 


k(s, t ) = cos(xs<) 

on [0,1]. A variety of parameters A that are close to the dominant characteristic values (the 
reciprocals of eigenvalues) are considered, as the equation becomes more difficult to solve as A 
approaches characteristic values. The dominant characteristic value that we use in our example is 
1.4278. The right hand function y(s) is so chosen that 

x(s) = e ;c cos(7s), 0 < s < 1 (46) 


Table I. The First Algorithm 


A 

Desired 

Estimated 

Actual 

Dimension (Level) 
Coarsest Finest 

1.00 

1.0E-6 

6.82E-7 

6.76E-7 

3(0) 

65 (5) 

1.40 

1.0E-4 

1 .62E-5 

1.60E-5 

5(1) 

65 (5) 

1.43 

1.0E-4 

1.31E-5 

1.31E-5 

5(1) 

129 (6) 
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In Table I, the Estimated column is computed using (39). As A approaches the characteristic 
value of 1.4278, both the coarsest grid level and the finest grid level were increased. In Table II, 
we give the iterative rate of convergence at each level, and the number of iterations performed at 
each level is also given in parentheses. As noted in section 3, only one iteration is needed as the 
level increases. Whenever only one iteration is performed at any given level, the iterative rate of 
convergence is the maximum contraction number £ in (19) estimated using (24) and (25). 


Table II. Iterative Rate of Convergence of The First Algorithm 


A 

Desired 


Level 


1 

2 

3 

1.00 

l.OE-6 

2.10E-2 (2) 

5.14E-2 (1) 

2.03E-3 (1) 

1.40 

1.0E-4 

2.10E-1 (2) 

5.3 IE- 2 (2) 

7.57E-3 (2) 

1.43 

1.0E-4 

- 

1.44E-1 (2) 

1.44E-2 (2) 



4 

5 

6 

1.00 

l.OE-6 

3.40E-3 (1) 

3.80E-3 (1) 

- 

1.40 

1.0E-4 

5.93E-2 (1) 

3.79E-3 (1) 

- 

1.43 

1.0E-4 

4.40E-2 (1) 

3.75E-3 (1) 

3.89E-3 (1) 


For the second algorithm, the coarsest level corresponds to two subintervals. In order to give a 
reasonable comparison with the first algorithm, we first give the results with 2 node points in each 
subinterval. Thus the quadrature order coincides with that of the first algorithm. 


Table III. The Second Algorithm with p=2 


Dimension (Level) 


A 

Desired 

Estimated 

Actual 

Coarsest 

Finest 

1.00 

l.OE-6 

6.82E-7 

6.76E-7 

4(0) 

64 (5) 

1.40 

1.0E-4 

1.62E-5 

1.60E-5 

4 (0) 

64 (5) 

1.43 

1.0E-5 

8.74E-6 

8.72E-6 

4 (0) 

128 (6) 


In the next table, we have results from the second algorithm with more node points on each 
subinterval. To show the superiority of the Gaussian quadrature rule, we give results for a smaller 
desired error for A = 1.40 and A = 1.43. 


Table IV. The Second Algorithm with p=3,4 
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A 

P 

Desired 

Estimated 

Actual 

Dimension (Level) 
Coarsest Finest 

1.40 

3 

1.0E-8 

1.95E-9 

1.93E-9 

6(0) 

96 (4) 

1.43 

3 

1.0E-8 

3.97E-10 

3.96E-10 

6(0) 

192 (5) 

1.43 

4 

1.0E-8 

6.52E-10 

6.28E-10 

8(0) 

64 (3) 


Table V. Iterative Rate of Convergence of The Second Algoritm with p-3, 4 


A 

P 

Desired 


Level 


1 

2 

3 

1.40 

3 

l.OE-8 

1.09E-4 (2) 

1.43E-6 (1) 

8.84E-4 (1) 

1.43 

3 

1.0E-8 

1.39E-3 (2) 

1.04E-2 (1) 

2.10E-4 (1) 

1.43 

4 

l.OE-8 

9.35E-6 (2) 

2.55E-3 (1) 

1.30E-5 (1) 





4 

5 

1.40 

3 

l.OE-8 


9.90E-4 (1) 

- 

1.43 

3 

l.OE-8 


2.36E-4 (1) 

2.42E-4 (1) 

1.43 

4 

l.OE-8 


- 

- 


CONCLUSION 


The piecewise Gaussian rule is superior to Simpson’s rule. However, as pointed out in section 
3, restrictions and prolongations are done with Nystrom type interpolation. And it involves 0(nf) 
multiplications at each level / without counting kernel evaluations. It appears that these 
operations cause the bottleneck of our algorithms. We are in the process of applying the idea 
suggested by Achi Brandt in [9] to our current algorithms which will reduce the operation count 
by far. Our preliminary results appear to be promising, and progress is being made in developing 

them further. 
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